Abelian Manna model on various lattices in one and two dimensions 
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We perform a high-accuracy moment analysis of the avalanche size, duration and area distribution 
of the Abelian Manna model on eight two-dimensional and four one-dimensional lattices. The results 
provide strong support to establish universality of exponents and moment ratios across different 
lattices and a good survey for the strength of corrections to scaling which are notorious in the 
Manna universality class. The results are compared against previous work done on Manna model, 
Oslo model and directed percolation. We also confirm hypothesis of various scaling relations. 
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1. INTRODUCTION 

Universality is the key feature of critical systems which 
justifies the analysis of (over-) simplified numerical mod- 
els of otherwise much more complex natural systems. 
When the concept of self-organised criticality (SOC) was 
introduced [1], it gained immediate popularity on the one 
hand because it attempted to explain the prevalence of 
scaling and fractality in nature and on the other hand 
suggested that many features of very different phenom- 
ena would be common to all of them, by the power of 
universality. At times disputed [2], it is now widely ac- 
cepted that the predictive power of SOC lies with its uni- 
versality; the elusive qualities of a critical point, normally 
confined to a very narrow range of a control parameter, 
are the norm in SOC models and shared between them. 

While universality has been at the centre of the debate 
about SOC [e.g. 3, 4], it comes as a great surprise that 
it has focused mainly on its occurrence across different 
models [e.g. 5], whereas very little work has been done 
on establishing it within the same models but across dif- 
ferent lattices [e.g. 6-9]. Recently, it has been shown [10] 
that on fractal lattices of same dimension but different 
structure, the Manna Model [11] behaves very differently 
The observation that critical phenomena and scaling dis- 
plays universal features on different lattices has tradition- 
ally been one of the most important insights, enabling, 
in particular, exact results [12, 13]. 

During the last decade or so, it has become increasingly 
clear that SOC is far more elusive than originally envis- 
aged. It is therefore all the more important to establish 
which features known from ordinary critical phenomena 



carry over to SOC. Independence of certain observables 
from the underlying lattice is one such aspect to substan- 
tiate. If universality was not to be found for the same 
model on different lattices, then the notion of universality 
in SOC as a whole had to be revised. 

In the following, the universality hypothesis in SOC 
is put to test by simulating the Abelian Manna model 
(AMM) [11, 14] on different one- and two-dimensional 
lattices. Owing to their robust scaling behaviour, in re- 
cent years, attention has focused on the AMM and the 
Oslo Model [15], which seem to be in the same universal- 
ity class [3]. This is contrasted by the poor scaling be- 
haviour of many other established SOC models, such as 
the Forest Fire Model [16, 17], the Bak-Tang-Wiesenfeld 
sandpile model [1, 18] or the Olami-Feder-Christensen 
Model [2, 19]. 

Apart from probing universality, corrections to scal- 
ing can also be tested, with the ultimate aim of find- 
ing a lattice that is particularly suited for simulating the 
AMM. In addition, rarely studied moment ratios are also 
reported which provide further support for the univer- 
sality of the AMM. This is of particular interest in one- 
dimension, where (logarithmic) corrections, are known to 
be remarkably strong [20]. 

The following is divided into four parts: In the next 
section, the model is defined and a few implementation 
details are discussed, including the method employed to 
determine the statistical error. All numerical results are 
collected in Sec. 3 and discussed in Sec. 4. 



MODEL DEFINITION AND METHODS 
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2.1. The Abelian Manna model on an arbitrary 
lattice 

The Abelian Manna model [11, 14, 21] is defined on a 
lattice C with N sites. Each site i on £ has neighbours, 
and is assigned to it a non-negative integer variable Zi 
which can be thought of as the number of particles at 
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that site or the height of a stack of particles there. The 
threshold height at all sites is 1 above which a site is 
said to be active or unstable, otherwise it is stable. The 
system evolves as follows. Driving: When the system 
is in a stable or "quiescent" configuration, i.e. Zi < 1 
for all sites, the system is charged by picking a site i 
at random and incrementing by 1. Relaxation: Ev- 
ery unstable site i relaxes by transferring two particles 
to its neighbours, possibly rendering the receiving site 
unstable. The recipient of each of the two particles is 
chosen randomly and independently. Unstable sites are 
updated in a random sequential order. By virtue of its 
Abelian nature the order of relaxations is irrelevant for 
(the statistics of) the final state of the lattice and the 
statistics of the avalanche size (see below). 

Each relaxation of a site constitutes a toppling, which 
in the bulk is conservative, i.e. the total J2i z i remains 
unchanged by bulk topplings. A (small) number of sites 
are considered dissipative boundary sites, which have 
"virtual" neighbours, that are included in the total count 
of neighbours ^ introduced above (but not in the num- 
ber of sites N). Virtual neighbours (or sites) never top- 
ple themselves, i.e. particles are lost from the system at 
receipt. Such sites therefore provide a dissipation mech- 
anism, and are introduced here solely as a bookkeeping 
device. As a general principle, the number of virtual 
neighbours of sites in the lattices discussed below are al- 
ways chosen so that the resulting at a boundary site 
matches that of a corresponding site in the bulk. In the 
lattices below, qi might take on at most two different 
values, yet bulk sites and boundary sites do not differ in 
that respect. 

An avalanche is the totality of all topplings until the 
system is quiescent again. The size s of the avalanche 
is measured as the number topplings performed between 
driving and quiescence, so that s = if no avalanche 
occurs after the driving. The area a of an avalanche 
is the number of distinct sites which received a particle 
during the avalanche. This includes the site charged by 
the driving, so that a > 1. The definition of time used to 
determine the duration t of an avalanche is based on the 
idea that each active site undergoes a Poissonian decay, 
i.e. all active sites topple with the same rate. This rate 
is chosen to be unity, so that on average time 1 /N a goes 
by until a site topples if N a sites are active. This is the 
amount by which time advances each time a site topples. 
Each time a site is picked at random, exactly two of its 
particles are redistributed. This procedure is, of course, 
not a faithful representation of a true Poisson process, 
which has random waiting times, but it is an increasingly 
good approximation in the limit of large activity N a [22] . 
Thanks to its Abelian nature [14], the model's evolution 
from quiescent state to quiescent state does not require 
a specific dynamics on the microscopic time scale, so the 
definition of time in the present case is to the same degree 
arbitrary as the dynamics in, say, model A [23] . 

The Abelian symmetry does not exist in the origi- 
nal model [11], but simplifies its implementation greatly. 



Originally, all particles were re-distributed at toppling 
which were subject to parallel updates, so that all sites 
that were active at time t had toppled before sites ac- 
tivated by a toppling during that time step. In that 
version of the model, it takes one time unit to update 
N a simultaneously active sites. In the Abelian version, 
this holds only for constant N a (number of active sites 
remaining unchanged). Yet, the update rate for both 
versions is l/N a when averaged over suitable intervals 
(one time unit in the former, one update in the later), 
although neither is Poissonian. Despite these differences 
in the definition, Dickman and Campelo [20] have shown 
that implementing parallel or sequential updating has no 
noticeable impact on the statistics. They also found that 
exponents derived from the Abelian variant of the Manna 
Model coincides with those published for the original ver- 
sion. 



2.2. The lattices 

In the following, we describe four one-dimensional 
(Fig. 1) and eight two-dimensional lattices (Figs. 2 and 
3). In the figures, links to adjacent sites are indicated by 
solid lines, those to virtual neighbours by dashed ones. 
As virtual neighbours are merely an accounting device, 
the positioning of the dashed lines in the figures is arbi- 
trarily chosen to match that of bulk links. 

The figures show in particular the boundaries, whose 
structure we maintained during finite size scaling, i.e. 
when we increased the lattice size, we made sure that 
the edges and corners of the larger lattice matched that 
of the smaller one. Initially we dismissed such issues 
as small corrections, but given the high accuracy with 
which we determine moments, changes in the boundaries 
became clearly visible, in particular in the Mitsubishi 
lattice (Fig. 3(d)). In hindsight this is hardly surprising 
given the great importance of transport of particles from 
the conservative bulk to the dissipative boundaries [24] . 

For some simple lattices, the number of sites N is ex- 
actly the power d (dimension) of their linear (Euclidean) 
extension L. In more complicated cases, this may hold 
only approximately. In order to maintain a reasonably 
high symmetry of the finite lattices, which have to be 
thought of as being cut out of an (infinite) tiling of the 
plane, we had to make compromises. In one dimension 
(d = 1, see Table I), A is a multiple of L for the first 
three lattices discussed. The forth, the Futatsubishi lat- 
tice (Fig. 1(d)), follows N = 3L + 1 for all sizes consid- 
ered. In two dimensions, d — 2, two different lengths L x 
and L y are introduced for the two different dimensions, 
defined in a way as natural as possible. The cutting of the 
finite two-dimensional lattices and the resulting shape of 
the boundaries was guided by the following principles, 
which clash and therefore hold only approximately: the 
number of sites in the lattice N — L x L y had to be as 
close as possible to the power of 2, the number of sites in 
a simple square lattice. The ratio L x /L y was to be held 
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TABLE I: Parameters for the one-dimensional lattices. Only the seven sizes eventually used are listed. The (nominal) linear 
size L, corresponding to the examples given in Fig. 1, is shown in brackets. 



lattice 


d 








N 








CPU 
time (h) 


Simple chain 1 


1024(1024) 


2048(2048) 


4096(4096) 


8192(8192) 


16384(16384) 


32768(32768) 


65536(65536) 


51350 


Rope ladder 


1 


2048(1024) 


4096(2048) 


8192(4096) 


16384(8192) 


32768(16384) 


65536(32768) 


131072(65536) 


54299 


nnn chain 


1 


2048(1024) 


4096(2048) 


8192(4096) 


16384(8192) 


32768(16384) 


65536(32768) 


131072(65536) 


7943 


Futatsubishi 


1 


3073(1024) 


6145(2048) 


12289(4096) 


24577(8192) 


49153(16384) 


98305(32768) 


196609(65536) 


54821 



TABLE II: Parameters for the two-dimensional lattices. Only the seven sizes eventually used are listed. The number of sites 
N equals the product L 



lattice 


d 




m x 


Cy 


m y 


L X X Ly 


CPU 
time (h) 


Square 


2 





1 





1 


256x256 


362x362 


512x512 


724x724 


1024x1024 


1448x1448 


2048x2048 


9530 


Jagged 


2 


1 


2 





1 


361x182 


511x257 


723x363 


1023x513 


1447x725 


2047x1025 


2895x1449 


9408 


Archimedes 


2 





4 





2 


360x182 


512x256 


724x362 


1024x512 


1448x724 


2048x1024 


2896x1448 


11106 


nc diagonal square 


2 


1 


2 


1 


2 


255x257 


361x363 


511x513 


723x725 


1023x1025 


1447x1449 


2047x2049 


9802 


Triangular 


2 





1 





1 


239x274 


337x389 


476x551 


673x779 


953x1100 


1347x1557 


1906x2201 


6305 


Kagome 


2 


1 


3 





1 


412x159 


583x225 


826x317 


1168x449 


1651x635 


2335x898 


3301x1271 


10881 


Honeycomb 


2 


1 


2 





2 


337x194 


475x276 


675x388 


953x550 


1347x778 


1903x1102 


2695x1556 


7842 


Mitsubishi 


2 


2 


3 


1 


2 


239x275 


338x387 


476x551 


674x777 


953x1101 


1349x1555 


1904x2203 


7573 



constant with increasing N. The (Euclidean) aspect ra- 
tio should be unity or very close, if that clashes with the 
shape of the boundaries. For each lattice m x , y and c XtV 
were set such that L XjV = c XiV (mod m Xiy ), to maintain 
the shape of the boundaries across different sizes. The 
parameters and resulting lattice sizes are listed in Ta- 
ble II. The number L y can normally be thought as a 
number of "layers" and L x as the number of sites in a 
layer. It is clear that such a definition, necessary because 
of the mismatch of the symmetry of lattice to the square 
symmetry of the cut-out, makes the linear sizes L x and 
Ly rather poor fitting parameters (see below). 

The size of each lattice is indicated exemplarily in the 
caption of its figure. The sizes of all lattices used in this 
article are listed in Tables I and II. In the following we 
describe the lattices in detail. 



2.2.1. One-dimensional lattices 

Simple chain (Fig. 1(a)): This lattice is the usual one- 
dimensional chain, where each site connects to two near- 
est neighbours. The leftmost and the rightmost sites have 
one virtual neighbour each. 

Rope ladder (Fig. 1(b)): The shape of this lattice is 
that of a rope ladder. It is a simple extension of the 
simple chain (Fig. 1(a)), which eventually leads to the 
square lattice (Fig. 2(a)). Each site has three nearest 
neighbours, one on the left, one on the right and one 
below or above it. Four boundary sites have one virtual 
neighbour each. 



Next nearest neighbour (nnn) chain (Fig. 1(c)): De- 
spite its triangular pattern, this is the simple chain 
(Fig. 1(a)) extended by allowing for next nearest neigh- 
bour interactions. This lattice was motivated by the ob- 
servation [25] that some models require such extensions 
in one dimensions to prevent degeneracy. Alternatively, 
it can be seen as the first step towards a two-dimensional 
triangular lattice (Fig. 3(a)). Each site has four neigh- 
bours, with boundary sites having either one or two vir- 
tual neighbours. 

Futatsubishi lattice (Fig. 1(d)): In appreciation of 
the Kagome lattice discussed below, "Futatsubishi" is a 
Japanese name, which translates to "two-diamond", re- 
flecting the shape of the lattice. The Futatsubishi lattice 
is fully contained in the Mitsubishi lattice (Fig. 3(d)) dis- 
cussed below (three of them meet at every vertex), or as 
a suitable slice of the square lattice (Fig. 2(a)) or the 
jagged lattice (Fig. 2(b)). On the present lattice, each 
site has either two or four neighbours in the bulk and 
two virtual neighbours at the boundary. 



2.2.2. Two-dimensional lattices 

Square lattice (Fig. 2(a)): This is the standard square 
lattice, the most commonly used lattice the AMM has 
been studied on. Each site has four nearest neighbours. 
As indicated by the dashed lines in Fig. 2(a), edge sites 
have one virtual neighbour and corner sites have two, 
leading to a constant qi — 4 for all sites i. 

Jagged lattice (Fig. 2(b)): This lattice is a square lat- 
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(a)The simple chain. 


L = 


W,N = 


10. 
































(b)The rope ladder, i 


I = 


10, N = 


20. 






\AAAAAAAAAZ" 



(c)The next nearest neighbour (nnn) chain. L = 10, N = 20. 



(d)Thc Futatsubishi lattice. L = 7, N = 22. 

FIG. 1: The four one-dimensional lattices considered in this 
article. Sites are shown as filled circled, adjacency is indi- 
cated by solid lines. Dashed lines indicate links to virtual 
neighbours. 



tice rotated by 45 degrees and fitted into a square shape 
which makes the boundary look jagged. The only differ- 
ence to the square lattice above is therefore the boundary, 
where sites have either two (edge) or three (corner) vir- 
tual neighbours, producing, again, a constant qi = 4 for 
all sites. There are various ways of cutting the lattice 
out of the bulk — we decided to maintain the left-right 
mirror symmetry, which results in three different bound- 
aries. In hindsight, a slightly different choice would have 
resulted in a lattice of higher symmetry and an aspect 
ratio closer to unity. The perfect match of the univer- 
sal features of the AMM on this lattice with those on all 
others is testament of the strong universal behaviour of 
the model. 

Archimedes lattice (Fig. 2(c)): The name of this lat- 
tice is normally complemented by (4, 8 2 ), which reflects 
the parameters for the general rule to construct these 
two-dimensional lattices [26]. In the present case, every 
vertex is surrounded by one square and two octagons. 
Each lattice site in the bulk has three neighbours. Along 
the boundary one of them is replaced by a virtual neigh- 
bour, so that qi = 3 throughout. 

Non-crossing (nc) diagonal square lattice (Fig. 2(d)): 
This lattice is based on a square lattice by adding alter- 
nate diagonals such that there are no crossing (nc) diag- 
onals. Each site has either four or eight nearest neigh- 
bours. Boundary sites have either one, three or five vir- 
tual neighbours. 

Triangular lattice (Fig. 3(a)): This lattice is the ap- 
proximate square-shaped clipping of a tessellation of the 
plane by triangles. It is probably the second most fre- 
quently studied lattice in statistical mechanics. Owing 
to its six-fold symmetry which clashes with the four-fold 
symmetry of the square, boundary conditions at the up- 
per and the lower edge differ from those on the left and 
on the right. This is not the case for three of the four pre- 
ceding two-dimensional lattices, but does similarly apply 
to the following three. In the present case, the number of 
virtual neighbours of sites along the edge varies between 




(a)The square lattice. (b)The jagged lattice. 

L x = L y = 6, N = 36. L x =A,L y = 9,N = 36. 




(c)Thc Archimedes lattice. (d)The non-crossing (nc) 
L x = 8, L y = 4, N = 32. diagonal square lattice. 

L x = L y = 5, N = 25. 

FIG. 2: The four two-dimensional lattices with four-fold sym- 
metry considered in this article. 



one and four, with a constant = 6 throughout. 

Kagome lattice (Fig. 3(b)): This lattice was first stud- 
ied by Syozi [12, 27]. Its name is Japanese, referring to 
the pattern of holes ("me") in a basket ("kago"). Each 
site in the bulk has four nearest neighbours. This lattice 
has a six-fold symmetry, which generates three different 
boundary conditions by the way we decided to cut it at 
top and bottom. Similar to the jagged lattice (Fig. 2(b)), 
in hindsight we may have picked slightly different bound- 
aries. 

Honeycomb lattice (Fig. 3(c)): Similar to the triangu- 
lar lattice, this lattice is a tiling by hexagons, leading to 
a honeycomb-shaped pattern. Each bulk site has three 
neighbours, = 3, and the number of virtual neighbours 
along the edge is one everywhere, even when top and bot- 
tom edges differ from those on the left and on the right. 

Mitsubishi lattice (Fig. 3(d)): This is Japanese which 
translates to "three-diamond" reflecting the shape of the 
lattice (the naming is inspired by the logo of the famous 
Japanese company of the same name). It is also known 
as "the diced lattice" [12, 27]. Each lattice site has ei- 
ther three or six nearest neighbours with a number of 
virtual neighbours at the boundary varying between one 
and four. Again, top and bottom edges differ from those 
left and right. 

Some of the two-dimensional lattices are related by 
transforms, which are frequently used in the analysis of 
equilibrium critical phenomena [12, 27-29]. Denoting the 
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(a)The triangular lattice. 
L x = 5, L y = 7, N = 35. 



(b)The Kagome lattice. 
L x = 10, L v = 4, N = 40. 





(c)The honeycomb lattice. 



(d)The Mitsubishi lattice. 
L x = 5,L y = 7,N = 35. 



FIG. 3: The four two-dimensional lattices with six-fold sym- 
metry considered in this article. 



duality transform by T>, the star-triangle transform by A 
and the decoration-iteration transform by I, the follow- 
ing equalities hold between two-dimensional lattices: 



2? (honeycomb) 
T> (Kagome) 
T> (Archimedes) 
2?(square) 
A (Mitsubishi) 
A(I(honeycomb)) 



triangular 

Mitsubishi 

nc diagonal square 

square 

triangular 

Kagome . 



in our experience produces very fast code with strongly 
optimising compilers. 

Storing the adjacency information in memory (rather 
than being implemented explicitly by rules), makes the 
code more flexible, but large lattices comparatively ex- 
pensive in terms of memory requirements. Significant 
amounts of memory memory are also required for the 
stack of active sites, which holds every site i whose Zi 
exceeds 1. Sites are placed on the stack at the time when 
they make the transition from Zi = 1 to z% = 2. Random 
sequential updating requires random access to that stack. 
Sites i picked from the stack topple only once and thus 
remain on the stack until z% < 1. At the time when an 
avalanche is triggered all sites might hold one particle, so 
that the theoretical maximum number of sites exceeding 
the threshold at any one time simultaneously is (7V+1) /2. 
Although this maximum is not reached in practice, build- 
ing in safeguards to protect smaller stacks from overflow- 
ing is computationally more costly than providing a stack 
as large as the theoretical maximum. 

A second stack is required to keep track of all sites top- 
pling during any one avalanche. Once a site is updated, 
a flag associated with it is changed, and the site's index 
is placed on the stack. It is an invariant that all sites 
with the flag raised are located on the stack. The area 
of an avalanche is the height of that stack at the end of 
the avalanche. The flags are reset by scanning through 
the stack. 

Even when using memory lavishly, memory require- 
ments for the AMM implementation as described above 
are rather modest compared to what any modern desktop 
computer has to offer. In one and two dimensions, large 
lattices are prohibitively large in terms of CPU time, not 
in terms of memory, as the average avalanche size grows 
quadratically in the linear system size, i.e. in the (aver- 
age) chemical distance to the open (dissipative) bound- 
ary. An additional, sometimes very noticeable constraint 
on systems with multiple cores or multiple logical cores 
(hyper-threading) is the memory bus, which can be alle- 
viated only partly by reducing the memory requirements. 



2.3. Implementation details 

Throughout this work we used the same implementa- 
tion in C of the AMM, which takes as an input the adja- 
cency information of the various lattices, that arc gener- 
ated separately from running the actual Manna Model. 
The full adjacency matrix is extremely sparse and it 
therefore makes little sense to store the information in 
that format. Rather, all sites are sequentially indexed 
and sites adjacent to a given site are listed by their in- 
dex in a sequence (of varying length). Negative indices 
indicate virtual neighbours. The adjacency information 
is then filled into a C struct for each site i, which con- 
tains a vector holding the sequence of indices of adjacent 
sites, their number, a flag whether the site has been hit 
by the currently running avalanche and finally the height 
Zi. Using indices rather than pointers to reference sites 



2.4. Output 

The output of the code described above is a string of 
moments, effectively subsamples, averaged over a num- 
ber of avalanches (ranging from many millions to several 
ten thousand), which we call "chunks" in the following. 
Typically, 100 to 10 000 chunks were generated for each 
lattice. The statistics of the chunks allows for an estimate 
of the statistical error, while a simple average (weighted 
by the chunk size if necessary) across chunks produces an 
unbiased, consistent estimate of the moments. 

The sizes of the chunks were chosen so that a new 
chunk would be produced every ten to sixty minutes. 
In one dimension the linear size of the lattices spanned 
about three orders of magnitude, in two dimensions 
the square root of that. The average size and thus 
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roughly the CPU time to produce a single avalanche 
grows quadratically in the linear size, ranging over six 
orders of magnitude in one dimension, and over three in 
two dimensions, see Tables I and II. The chunk sizes 
have to be adjusted accordingly. At the same time, 
the avalanche- avalanche correlation time increases like a 
power of the linear extension, but not with the dynamical 
exponents z, which defines the link between microscopic 
and macroscopic time scale, but with L D ~ d , which is a 
measure of the characteristic fluctuations of the AMM in 
the interface picture [30-32]. 

The chunks were generated on the SCAN facility at 
Imperial College London, which harvests CPU time from 
undergraduate computing facilities when not used by stu- 
dents, providing up to 706 logical CPUs simultaneously, 
mostly in the form of Intel® Core2 processors with 
2.66 GHz. In this setup, all I/O is done over the net- 
work, which in principle constraints the amount of out- 
put per process running. Given the drastic difference in 
CPU time required to generate a single avalanche, by far 
the most CPUs were assigned to the largest systems. In 
two dimensions, the smallest four system sizes each had 
one logical CPU to themselves. As a result, the small- 
est systems (that we kept, see below) have sample sizes 
exceeding 10 9 avalanches, which leads to highly accurate 
estimates for their momenta. In total, about 168 415 and 
72 451 hours of CPU time were spent on generating the 
data in one and two dimension, respectively. 

For simplicity, checkpointing was implemented only 
right after a chunk was written out. The power cycling of 
the computer setup thus limited the amount of time avail- 
able to generate a single chunk. While the original inten- 
tion was to choose the chunk size such that correlations 
(which normally come in the form of anti-correlations) 
are negligibly small, this turned out to be unsustainable 
for the very large system sizes. However, uncorrelated 
chunks greatly facilitate the calculation of statistical er- 
rors, compared to, say, a full-blown resampling plan [33]. 
To this end, chunks were merged during post-processing, 
as discussed below. 

As each instance of the AMM was started with an 
empty lattice, a generous amount of chunks was dropped 
as transient from the set considered in the subsequent 
data analysis. This equilibration "time" was estimated 
by inspection of individual series of chunks, as a multi- 
ple of the time to "obvious" stationary. For the largest 
lattices in one dimension typically 10 5 avalanches (more 
for larger lattices) were rejected as transient and 10 6 re- 
tained for statistics in one dimension. In two dimensions, 
typical numbers were 6 • 10 6 as transient and 400 • 10 6 for 
statistics [69]. As a rule of thumb, at least 3/2 times the 
number of sites in the lattice avalanches were removed as 
transient. For smaller lattices a much larger fraction was 
taken, as for them, equilibration was often apparently 
reached within a single chunk. Equilibration can also be 
observed in the density of particles, which, at least in 
one dimension, initially grows almost perfectly linearly 
with only a minute amount of dissipation. Once a cer- 



tain fraction of sites is occupied, the occupation density 
displays very litte relative fluctuations, for large lattices 
of around 2 • 10~ 4 . 

The transient serves the additional purpose of warm- 
ing up the random number generator (RNG), which was 
of course seeded uniquely in each instance, except for the 
Mitsubishi lattice, whose 190 different seeds equal those 
used for honeycomb lattice. All results presented in the 
following are based on the Mersenne twister [34], which 
has received some criticism for its correlations across dif- 
ferently seeded instances [35]. The independence of the 
present results from the RNG was tested by re-running 
a few setups with Marsaglia's KISS RNG [36]. 



2.5. Statistical error 

All results presented below are based on moments of 
the avalanche size, area and duration. 

The first moments of all three observables generally 
posed a problem, possibly because they were determined 
with an accuracy so high that it was virtually impossible 
to account for their corrections to scaling (see below). 
The statistical error of all moments was calculated on the 
basis of chunks whose size was chosen as to ensure their 
independence. This was tested by firstly calculating their 
autocorrelation function and, secondly, by successively 
increasing their size, i.e. merging them, to probe whether 
the statistical error derived on their basis was affected by 
this operation. In fact, for the largest system sizes (and 
thus smallest chunk sizes) tested, some correlations were 
visible (correlation length of about 0.7 chunks), and we 
decided to merge ten consecutive chunks throughout. It 
turned out, however, that only the statistical errors of 
the first moments were noticeably, yet still insignificantly, 
affected at all by this operation. 

It is straight-forward to determine estimates and their 
estimated variance on the basis of independent chunks. 
Denoting the observable in the zth chunk by c,, an unbi- 
ased, consistent estimator [37] of its population average 
(c) is 



— y 



(1) 



given a sample size of M c . In the following, we will de- 
note that estimate itself by (c). An unbiased, consistent 
estimator of the variance of this estimator is 



1 

Ml 




(2) 



These estimates were used as input values for the estima- 
tion of the exponents, based on the scaling assumption 
discussed in the following. For completeness, covariances 
of the estimators for the averages of two observables c 
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and c', as used in Eq. (18), were estimated using 



cov(c, c') 



M„- 1 



1=1 



i=l 



i=l 



(3) 



As detailed below, the finite size scaling of the moments 
of each observable was determined in the form 



,(*) 



(x n ) = a^„L M " + corrections 



(4) 



where x stands for the size s, the area a or the duration 
t. The combined information of the scaling of a num- 
ber of different moments n was then used to estimate 
the finite size scaling exponents. Even when we imple- 
mented the first eight moments, we used only moments 
2 to 4 {e.g. 

/•4 S ^> r 1 ^) when fitting the avalanche 

dimension (using 2 = (2 — t)D to determine r) and mo- 
ments 2 to 5 to fit exponents characterising the scaling 
of avalanche area and duration (all exponents are de- 
fined below). Different moments estimated on the basis 
of a single Monte-Carlo simulation are not independent. 
To account for that, we (rather generously) multiplied 
the statistical error of each moment by the square root 
of the number of moments considered simultaneously, as 
if each moment was determined independently from the 
others [38] . For example, the statistical error of the sec- 
ond, third and fourth moment of the avalanche size was 
multiplied by \/3, before fJ,^\ and jif^ respectively 
were determined. It seems that these correlations are 
frequently ignored in the literature. This procedure does 
not account for correlations in the estimates of moments 
of different observables and thus our results for the fi- 
nite size scaling exponents for different observables are 
not independent. Multiplying again by the square root 
of the number of different observables considered, may, 
however, seriously overestimate the impact of these cor- 
relations. 



2.6. Finite size scaling 

Making the usual finite size scaling assumption for ob- 
servable x, its probability density function (x) fol- 
lows 



V (x) (x;L) = a x x- 



b x L D * 



(5) 



asymptotically in large x ^ x with lower cutoff xq, lin- 
ear system size L, non- universal metric factors a x and b x 
and universal exponents t x and D x . Below, the fractal 
dimension of avalanches is denoted by D a [39] and the 
avalanche area exponent by r a . For historic reasons, the 
dynamical exponent is denoted by z (rather than D t ), 
the avalanche duration exponent by a (rather than r t ), 
and finally the avalanche dimension by D (rather than 



D s ) and the avalanche size exponent by r (rather than 
t s ). The universal, dimcnsionless scaling function Q x of 
a dimensionless argument decays, for large arguments, 
faster than any power law, so that all moments 



rOO 

x n ) (L) = / dxV (x) {x;L)x n 
Jo 



(6) 



exist for any finite system and n > 0. Provided that 
n + 1 — r > one can easily show [40, 41] that gap 
scaling follows, so that 

POO 

(x n )(L) = a x (b x L D °) n+1 - Tx dyy n -^g x (y) (7) 

Jo 

asymptotically in large L, or, more accurately, 



< km - J; / v ' < oo . 



(8) 



To determine D x (and r x , if independent), the leading 
order scaling of the moments according to Eq. (4) was 

fx) 

estimated, fitting the resulting exponents fin against 
D x (n+ 1 —t x ), without allowing for any further correc- 
tions. The results are shown in Tables III and IV. 

Because of particle conservation in the bulk, every 
particle placed on the lattice anywhere can only leave 
through the boundary. On the way there, it performs 
an independent random walk, even when it occasion- 
ally rests [3]. Every move by a particle is caused by a 
site's toppling and the number of particle moves during 
an avalanche is therefore exactly twice the number of 
topplings. At stationarity one particle leaves the system 
per particle added and (half) its trajectory length is its 
contribution to the various avalanches it has been part 
of. The average contribution per particle and thus per 
avalanche is the average trajectory length, which scales 
like L 2 independent of the dimensionality of the lattice 
[13] and many of the details of the boundary. This ar- 
gument remains valid if avalanches of size are excluded 
from the average (which we did not), provided the prob- 
ability of producing an avalanche of size (i.e. hitting an 
empty site) does not converge to 1. 

As a result (s) cx L 2 asymptotically, or, here, (s) cx 
N 2 / d , i.e. fj,^ = 2 and under the assumption of gap 
scaling 2 = (2 — t)D. This identity has been used in 

(s) 

the fitting of the avalanche dimension, i.e. /i„ for n = 
2, 3, 4 was fitted against 2 + D(n — 1). At the same time, 
comparing the estimate for (Tables III and IV) to the 
exact value 2 allowed us to assess the fitting procedure, 
in particular the form of the corrections discussed below. 
Defining 



-Ex = D x (t x - 1) 



(9) 



the assumption of a sufficiently [42, 43] narrow joint prob- 
ability density function [39, 44, 45] produces the scaling 
law S := E s = E t = E OJ i.e. 



-E = D(t - 1) = z(a - 1) = D a (r a - 1) . 



(10) 
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As this is not a mathematical identity [30] (in fact, it 
seems to be broken in the original Manna Model [39]) it 
is tested as a scaling hypothesis below. The exponent E 
can be seen as a replacement for the exponents t, a and 
T a , as [in^ = nD + E, ^ = nz + Y<, ^ = nD a + E for 
n> t — 1, n> a — 1 and n > r a — 1 respectively. 

Finally, as in the BTW model [46] , in the Manna Model 
avalanches are often assumed to be compact [4, 42], i.e. 
D a = d (supposedly up to the upper critical dimension, 
where we expect dangerous irrelevant fields to spoil this 
scaling relation) . There is no mathematical proof for this 
feature, yet numerically it is well verified, see Tables III 
and IV. Clearly D a < d, which served as another con- 
straint to assess the quality of the exponents extracted. 



2.6.1. Corrections to scaling 
(x) 

The exponents /x„ characterise the asymptotic scal- 
ing of the moments in large L. It is widely known [47], 
however, finite lattices suffer from finite size corrections, 
which manifest themselves as sub-leading terms to be in- 
cluded on the right hand side of Eq. (4) [42]. A priori, 
the structure of the corrections is not known, yet they 
have a marked impact on the quality of the results, as 
they are imposed when fitting the data. 

Given that there is often no natural way of defining 
the linear extent of a lattice, we decided to replace L (as 
in Eq. (4)) by N x / d . We considered a host of different 
fitting functions, such as 



(x n ) (N) = A x , n N^'' d (11) 

(x n ) (N) = A x ^N^ x) /d + Bx nN ^/d-^, n /d {12) 

(x n )(N) = A^N^^ + B^N^/d-m (13) 

(x n )(N) = A^N^^ + Bx^N^'*- 1 ' 2 (14) 



and eventually settled for 



(x n ) (N) = A x , n N^ ] / d + B^N^/d-^ 

+ ^„A^-i/2 (15) 

which yielded particularly good estimates. In particular, 

(s) 

fi\ — 2 was reproduced quite reliably. The quality of 
the estimates was assessed on the basis of the goodness of 
fit determined in the Levenberg-Marquardt least square 
fitting routine [48] . For the vast majority of moments and 
lattices, we could have dropped the last term in Eq. (15) 
arriving Eq. (13) and still achieved a goodness of fit (q- 
value) of more than 0.9. However, the first moments of 
the avalanche size, whose finite size scaling exponent is 
the only exactly known one, was particularly poorly fitted 
without that term. For consistency, we decided to fit all 
moments using Eq. (15), achieving typically q- values of 
better than 0.9, suggesting that we overestimated the 



statistical errors. In all result tables, fits that had a q- 
value of less than 0.1 are marked as such. 

To reduce the impact of a possible dependence on the 
(arbitrary) choice of the initial values of the free param- 
eters, the number of terms in the fitting function was 
increased successively starting from Eq. (11), using the 
estimates of the parameters in the previous fit as the ini- 
tial value of the same parameters in the next. 

The simplicity of Eq. (15) meant that we had to drop 
results for small system sizes, which suffer from stronger 
finite size corrections, yet were determined with much 
greater accuracy than those of the bigger systems. This 
is a common theme in the present work: Moments in 
small system sizes were determined with such great accu- 
racy that very many (a priori unknown) correction terms 
would have to be included to account for all such details. 
At the same time, it makes little sense to have almost as 
many free parameters in the fitting function as there are 
data points to fit. In order to retain the goodness of fit 
as a meaningful device to determine the quality of the fit, 
we therefore removed the smallest four system sizes from 
the procedure (in one and two dimensions), keeping the 
seven system sizes listed in Tables I and II. 

Increasing the system size in order to suppress correc- 
tion terms comes at the price of increased relative error 
if ^ > 1. According to Eq. (4) and Eq. (8), the vari- 
ance of the nth moment has leading order L r>x ^ 2n+1 ~ Tx \ 
i.e. the relative error scales like L Dxt ^ Tx ^ 1 ^ 2 . Moreover, 
correlations are expected to die off after L D ~ d 7 which 
reduces the number of effectively independent measure- 
ments with increasing system size. 



2.7. Moment ratios 

Ratios of products of moments, which (to leading or- 
der) are independent from L, characterise the scaling 
function Q x , Eq. (5), directly. In equilibrium phase tran- 
sition, the so-called Binder-cumulant [49, 50] is the best 
known such ratio, signalling the deviation from a Gaus- 
sian distribution of the order parameter around the crit- 
ical point. There are many ways of constructing suitable 
moment ratios; assuming (x n ) cx L D *( n + 1 - T *) ; it is easy 
to see that any ratio of products of moments, which has 
the same number of moments (to cancel t x ^ 1) and 
the same sum of orders of moments (to cancel D x ) in 
numerator and denominator leads to a non-scaling quan- 
tity. Since the second moment is positive and bounded 
away from 0, traditionally moment ratios are formed by 
dividing by a power of it. Moreover, in many phase tran- 
sitions, the order parameter follows a distribution with 
t x = 1, which removes the constraint of having the same 
number of moments in the numerator and the denomina- 
tor. 

While the sets (x n - ,n ) (x n+m ) / (x n ) 2 attracts by its 
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simplicity and symmetry, the set 



„(x) = (*> 



n-2 



<X 2 > 



2\n-l 



(16) 
} = 1 



has the particularly nice feature that g{ ' — g' 2 
by definition, which fixes the metric factors a x and b x in 
Eq. (5) by imposing for n = 1,2 

/■OO 

ff <*> = / dyy n -^g x (y) , (17) 
Jo 



see Eq. (7) which is then consistent with Eq. (16) for all 
n. 



The statistical error a (dr?^ °f the estimator of 
Eq. (16) is to leading order in the sample size given by 



, COv(x , X ) 

«*» 2 



+ (n-l) : 



, cov(x 2 ,x 2 ) 



«x 2 » 



_,cov(x 2 , x 1 ) . _cov(x 2 ,x") n/ _ cov(x 1 ,x") . 
2(n - 2)(n - 1)-^-^ - 2(n - 1) , 2 \ + 2(n - 2) ^ ^ ] , (18) 



(x) (x 2 > 



(x 2 ) (x n ) 



(x) (x«) 



which matches perfectly (typically the first three or four 
significant digits) the error as found by the subsampling 
using chunks, i.e. determining g^ for each chunk and 
estimating the error by the square root of its variance 
over the number of chunks. In Eq. (18) (x™) strictly de- 
notes the estimator of the nth moment and cov(x™,x m ) 
the estimated covariance of the nth and mth moment, 
see Eq. (3). Consistent with the preceding discussion, we 
used averages and statistical errors derived from chunks. 



All lattices were set up with the intention of creating an 
aspect ratio of 1, which is trivial as long as the lattice has 
a four-fold symmetry. In particular for lattices without 
that symmetry, such as the triangular lattice (Fig. 3(a)), 
the Kagome lattice (Fig. 3(b)), the honeycomb lattice 
(Fig. 3(c)) and the Mitsubishi lattice (Fig. 3(d)), but 
also, say, the jagged lattice (Fig. 2(b)), the aspect ratio 
might deviate slightly from unity and converge to 1 only 
with increasing system size. In any case, the aspect ratio 
might be more reasonably be defined using the Manhat- 
tan distance across the lattice. It is well known that uni- 
versal scaling exponents are generally independent from 
the aspect ratio, whereas finite size scaling functions are 
not [51]. Therefore, deviations of the moment ratios in 
particular in case of the lattices listed above are expected 
(but did in fact not materialise). Surprisingly, even when 
there is every reason to assume that no such problem can 
occur in one dimension, their moment ratios proved par- 
ticularly difficult to analyse. 



3. RESULTS 

3.1. Avalanche exponents 

After stripping off the transient, merging chunks as 
discussed above, deriving average moments and errors 
using the procedures described above, the scaling expo- 
nents (iffi where fitted using Eq. (15) for the three dif- 
ferent observables, x = s,t,a, avalanche size, duration 
and area respectively. As mentioned above, correlations 
between moments were taken into account by multiply- 
ing the statistical error of the moment by the square root 
of the number of moments considered. Allowing for no 

further corrections, [in^ were fitted for each of the twelve 
different lattices separately against 

/£> = 2 + D(n - 1) (19) 

for n = 2, 3, 4. The results are collected in Table III. The 
scaling law 2 = (2 — t)D used in Eq. (19) was probed 

independently; fi[ s ^ = 2 is a mathematical identity, but 

(s) 

valid only asymptotically and the deviation of fi\ from 2 
can therefore serve as an indicator to assess the quality of 
the fitting routines and the data and can help confirming 
that "asymptotia is reached". The estimate for ^ on 
the basis of Eq. (15) is shown in Table III alongside the 
other finite size scaling exponents. The exponent r stated 
in Table III is derived from the estimate of D using r = 
2-2/D. 

For all observables (size, area and duration), the first 
moments turned out to be problematic. According to 
Eq. (6) moments n < r — 1 remain finite in the ther- 
modynamic limit, consistent with our observation that 
smaller moments generally require more correction terms. 
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The average avalanche size is particularly difficult to han- 
dle, which is determined, due to the presence of anti- 
correlations [31, 32], with incredible precision, typically 
with a relative error of the order 10 ~ 5 . We therefore de- 
cided to omit first moments from the determination of 
the finite size scaling exponents throughout. 

For x = t and x — a no scaling laws were used (even 
when D a = d is generally assumed to hold) and so the 
exponents fi^ were fitted against: 

- z(n + l-a) (20) 
fi^ = D a (n + l-r a ) (21) 

for n = 2,3,4,5. 

Fitting moments beyond n = 5 proved very difficult. 
We decided to drop all moments beyond the fourth for 
avalanche size and beyond the fifth for the avalanche area 

(x) 

and duration, as fi n ; became clearly dependent on the 
choice of the initial values of the fitting parameters. As 
mentioned above, in all fitting schemes used, we increased 
the number of free parameters successively and used the 
estimates of the fitting parameters of one scheme as the 
initial values for its extension. For example, we used A x>n 

and fin^ /d from a fit against Eq. (11) as initial values in 
a fit against Eq. (13), which in turn produced the initial 

values of A x>n , £? X; „ and fifi^ jd to fit with Eq. (15). We 
observed this procedure in all fitting schemes discussed 
below. 

As the variance of the nth moment scales like 
_L_D x (2n+i-r x \ its numerical estimate is increasingly af- 
fected by the floating point precision (double-extended 
throughout) — equivalently, the typical largest measure- 
ment of the nth moment scales like the nth power of the 
cutoff, L nD ^ which for n = 6, D x = 2.25 and L = 2 16 
is 2 216 . Given that the smallest event size is 0, this is 
to be compared to the 64 bits in the mantissa on a long 
double on the x86 architecture. 

The exponent T, x can be derived either from the def- 
inition Eq. (9) through D x and t x or by independently 

(x) 

fitting fi n ; against nD x + . We took that approach for 
x = t and x = a. Using the three different observables 
size, duration and area, provides effectively estimates for 
D(t — 1), z(a — 1) and D a (r a — 1) respectively. In the 
case of D(t — 1) this is in fact exactly D — 2, since we 
imposed D{2 — t) = 2 when estimating D. The entry for 
E s in Table III is therefore derived from the estimate for 
D. Except for r, all other entries in Table III are based 

(x) 

on fitting fi n directly. 

Table III provides very strong support for universal- 
ity across different lattices. Under the assumption that 
universality holds, estimates for exponents gained from 
different lattices can be taken together to produce an 
overall estimate. The result of that procedure is shown 
in Table IV. As discussed further below, we regard only 
the data referring to the fitting function Eq. (15), shown 
in bold, as reliable, producing the most consistent and 
robust results. 



Fits with Eq. (15) produced very large g-factors, sug- 
gesting we had been too generous either with estimating 
the statistical errors or with the number of free param- 
eters. We therefore also tried Eq. (13), which, however, 
gave partly inconsistent results. In particular estimates 
for fi[ s ^ deviated from the exact value 2 by about 30 
standard deviations (although the relative error was only 
3T0~ 3 ). For the simple chain the avalanche size moments 
tested displayed a poor quality of fit using Eq. (13), as 
did the first moments of the avalanche size for all lattices 
except for the nnn chain (Fig. 1(c)), the Futatsubishi 
lattice (Fig. 1(d)) and the Archimedes lattice (Fig. 2(c)). 
The results for the fits against Eq. (13) are also sum- 
marised in Table IV. A bracket [•] indicates finite size 
scaling exponents not being fitted with a goodness of fit 
better than 0.1. 

In general, the Futatsubishi lattice and the simple 
chain were particularly difficult to fit using whichever 
fitting function. 

In order to extract the exponent of the sub-leading 
terms, the remainder (x n ) (N) — Au^iV"'^ was fitted 
against 

B Xtn N^ x) / d -^,n/d (22) 

determining 0J x>n . This procedure was aiming much 
more at a qualitative result rather than a quantitative 
one and generated rather noisy estimates. The Futat- 
subishi and the triangular lattices proved particularly 
difficult to handle. The data in d = 1 produced, un- 
fortunately quite inconsistently w s ,n ~ 0.28, cot.n ~ 0.20 
and Lo atn rts 0.20, while d = 2 produced, slightly more 
consistently uj St . n w 0.23, uj t , n ~ 0.32 and u> a . n s=s 0.47 
fairly independent of lattice and n (but more reliably for 
large n and observables other than the avalanche size). 
These exponents could in turn be used in Eq. (12) to fit 

(x) 

the data for fi n at fixed u} x>n . The resulting overall esti- 
mates for the finite size scaling exponents are also shown 
in Table IV. 



3.2. Moment ratios 

Similar to the plain moments, one has to allow for 
corrections when fitting moment ratios. Most two- 
dimensional lattices (except nc diagonal square lattice 
(Fig. 2(d)) and Mitsubishi lattice (Fig. 3(d))) produced 
consistent results with a goodness of fit of greater than 
0.1 with 

g { : ] + D Xin N~ - 25 (23) 

but in order to capture all lattices and for consistency 
with the above we decided to add a further correction, 
finally fitting against 

gW + D x , n N~ - 25 + E x , n N-°- 5 . (24) 

(x) 

In contrast to the finite size scaling exponents fi n of 
the moments considered above, all moment ratios were 
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TABLE III: Summary of all exponents characterising the avalanching in the twelve different lattices, using Eq. (15). 
estimates for r and D(r — 1) were not determined by fitting the data, but through the scaling relation D(2 — t) = 2. 



The 
The 



estimates for /ij verify this scaling relation. The estimates in the last three columns should coincide under the narrow-joint- 
distribution assumption, Eq. (10). Estimates for different observables are not independent. 

lattice d D t z a D a 



0) 
Mi 



Simple chain 
Rope ladder 
nnn chain 
Futatsubishi 



1 2.27(2) 

1 2.24(2) 

1 2.33(11) 

1 2.24(3) 



1.117(8) 
1.108(9) 
1.14(4) 
1.105(14) 



1.450(12) 
1.44(2) 
1.48(11) 
1.43(3) 



1.19(2) 
1.18(3) 
1.22(14) 
1.16(6) 



0.998(4) 
0.998(7) 
0.997(15) 
0.999(15) 



1.260(13) 
1.26(2) 
1.27(5) 
1.24(5) 



2.000(4) 
1.989(5) 
1.991(11) 
2.008(11) 



0.27(2) 
0.24(2) 
0.33(11) 
0.24(3) 



0.27(3) 
0.26(5) 
0.3(2) 
0.23(9) 



0.259(14) 
0.26(2) 
0.27(5) 
0.24(5) 



Square 

Jagged 

Archimedes 

nc diagonal square 

Triangular 

Kagome 

Honeycomb 

Mitsubishi 



2 2.748(13) 
2 2.764(15) 
2 2.76(2) 
2 2.750(14) 
2 2.76(2) 
2 2.741(13) 
2 2.73(2) 
2 2.75(2) 



1.272(3) 
1.276(4) 
1.275(6) 
1.273(4) 
1.275(5) 
1.270(4) 
1.268(6) 
1.273(6) 



1.52(2) 
1.54(2) 
1.54(3) 
1.53(2) 
1.51(2) 
1.53(2) 
1.55(4) 
1.54(3) 



1.48(2) 
1.49(3) 
1.50(3) 
1.49(2) 
1.47(3) 
1.49(2) 
1.51(4) 
1.50(4) 



1.992(8) 

1.995(7) 

1.997(10) 

1.992(7) 

2.003(11) 

1.993(8) 

1.990(13) 

1.999(12) 



1.380(8) 

1.384(8) 

1.382(11) 

1.381(8) 

1.388(12) 

1.381(9) 

1.376(14) 

1.387(12) 



1.9975(11) 

2.0007(12) 

2.001(2) 

2.0005(12) 

1.997(2) 

1.9994(12) 

2.000(2) 

1.998(2) 



0.748(13) 

0.764(15) 

0.76(2) 

0.750(14) 

0.76(2) 

0.741(13) 

0.73(2) 

0.75(2) 



0.73(4) 
0.76(5) 
0.78(6) 
0.75(4) 
0.71(6) 
0.75(5) 
0.79(8) 
0.77(7) 



0.76(2 
0.77(2 
0.76(3 
0.76(2 
0.78(3 
0.76(2 
0.75(3 
0.77(3 



TABLE IV: Overall estimates of scaling exponents in one and two dimensions. Only the fits using Eq. (15), based on the data 
in Table III and shown in bold, are fully reliable. Entries for Eq. (13) and Eq. (12) are for comparison to other estimates only. 
Fits with a goodness of less than 0.1 are marked by [•]. The estimate for E, Eq. (10), is based on all estimates for D(t — 1), 
z(a — 1) and D a {r a — 1) in Table III. Their correlation is taken into account by multiplying their respective error by \/3- 

d function D r z a D a r a — £ 

1 Eq. (15) 2.253(14) 1.112(6) 1.445(10) 1.18(2) 0.998(3) 1.259(11) 0.26(2) 

1 Eq. (13) [2.265(4)] [1.117(2)] [1.449(2)] 1.172(3) 1.0000(6) 1.249(2) 0.249(3) 

1 Eq. (12) [2.2520(3)] [1.11188(H)] [1.4632(6)] [1.219(2)] 1.0000(8) 1.276(2) [0.297(3)] 

2 Eq. (15) 2.750(6) 1.273(2) 1.532(8) 1.4896(96) 1.995(3) 1.382(3) 0.761(13) 

2 Eq. (13) 2.7698(12) 1.2779(3) 1.5407(14) 1.498(2) 1.9990(5) 1.3843(6) 0.768(2) 

2 Eq. (12) [2.7673(3)] [1.27728(7)] 1.541(2) 1.501(2) 1.9985(6) 1.3853(6) 0.770(2) 



fitted as if they were independent, i.e. considering them 
simultaneously may be misleading as they are correlated 
and these correlations have not been accounted for. 

The results are shown in Table V. In general, two 
dimensional lattices are much better than one dimen- 
sional ones. The observable most easily fitted is the area 
size distribution (which might be caused by the plain 
amplitudes A a ^ n being universal, see below). The only 
two-dimensional lattice that displayed low goodness of 
fit throughout was the triangular lattice, while the hon- 
eycomb lattice had a single poorly fitting ratio. In one 
dimension, the picture is reversed; There is hardly any 
reasonably fittable moment ratio. Fits which produced 
a goodness of less than 0.1 are marked in Table V again 
by [•]. Results are rather noisy for the highest moment 
ratios, which might suggest an explanation for the slight 
inconsistencies, which are not covered by the statistical 
error, for example for in nc diagonal square lattice 
(Fig. 2(d)) and Mitsubishi lattice (Fig. 3(d)). Yet, in one 
dimension, it is the lower order moment ratios that were 
most difficult to handle. 

The overall estimates for the moment ratios are shown 
in Table VI. We refrained from stating an overall esti- 



mate, where not at least two lattices produced reliable 
estimates as shown in Table V. Those that show low 
quality of fit are marked as above. 

There is, again, a certain sensitivity to the fitting func- 
tion; Eq. (23) gives slightly incompatible results, with 
remarkably small error bars. Given what has been said 
about the goodness of fit, we place our confidence in the 
results presented in Table VI. 



4. DISCUSSION 

Exponents (Table III) and, at least in two dimensions, 
moment ratios (Table V) are universal across different 
lattices. Subleading orders of moments are noisy, but 
still fairly consistent. There can be little doubt that the 
Abelian Manna model displays all the hallmarks of a crit- 
ical system as they are known from equilibrium critical 
phenomena. The exponents shown in bold in Table IV 
and the moment ratios in Table VI (except those shown 
in brackets) are perfectly consistent across all our results 
and represent reliable, high accuracy estimates of the uni- 
versal quantities characterising this universality class. 
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TABLE V: Estimates for the moment ratios as defined in Eq. (16) obtained by fitting the relevant ratios against Eq. (24). Fits 
with a goodness of less than 0.1 are marked by [•]. 
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In the literature can be found a large number of es- 
timates for the exponents in Table IV. Traditionally, 
the AMM is studied in two dimensions and therefore 
many more estimates are available in that dimension. 
Given that the AMM is in the same universality class 
as the Oslo Model [3, 15], there is a second source for 
comparison. What makes the comparison more compli- 
cated is the fact that many authors have studied variants 
of the Manna Model (in fact, the Abelian version stud- 
ied here is a variant of the original model); for example 
Dickman and Campelo studied the Manna Model with 



height restrictions [20] and Liibeck and Heger studied 
its "fixed energy sandpile" version [56, 60]. Table VII 
collects a broad range of estimates across the literature, 
which nevertheless provide a perfectly consistent picture. 
The present work clearly improves on comprehensiveness 
and on accuracy, which for most estimates is improved 
by one digit. The fact that some estimates in the liter- 
ature have even smaller error bars than ours might be 
partly due to our over-estimation of statistical errors but 
also due to other authors using models that are better 
behaved in one dimension, as in [52, 53]. 
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TABLE VI: Overall estimates for the moment ratios defined 
in Eq. (16), based on the data presented in Table V, if enough 
data was available. 



1 s - 
1 t - 

1 a [1.3320(5)] 1.961(2) 3.035(4) 4.838(9) 

2 s 1.8273(14) 4.363(10) 12.32(6) [39.3(3)] 
2 t 2.11423(98) 5.939(8) 19.88(6) 75.8(4) 
2 a 1.7501(4) 3.711(2) 8.699(8) 21.72(3) 

Table VII also contains the avalanche exponents for the 
directed percolation (DP) universality class. The expo- 
nent D is derived through the identity D = d + z — (3/v±, 
where —(3/v± is the finite size scaling exponent of the 
activity in DP [70]. The dynamical exponent z is some- 
times used in the DP literature as what would have been 
in the present notation 2/z. The exponent for D a is 
based on the asumption of compact avalanches [55]. Al- 
though sometimes disputed in the literature [61, 62] using 
avalanche exponents leaves little doubt that the Manna 
universality class differs from DP. As DP is normally 
performed with periodic boundaries and with different 
observablcs, there is, to our knowledge unfortunately 
no published work on the moment ratios we considered 
here. In fact, depending on the definition of the ensemble 
[56, 63-65], some moment ratios in DP can be undefined. 

Certain finite size scaling features specific to SOC are 
confirmed as well: Compactness of avalanches, D a = d 
is very strongly supported (in two dimensions the devi- 
ation of D a from d is slightly bigger than one standard 
deviation, though), as is the universality of S, Eq. (9). 
It is reassuring that the asymptotics of the first moment 

(s) 

of the avalanche size, = 2, are recovered, validating 
our numerical schemes. 

Although we invested more than twice as much CPU 
time (in absolute terms) in one dimensional lattices, the 
results are significantly noisier than for two dimensional 
lattices. Providing a sufficient number of correction 
terms still produces consistent data (Table IV), but the 
results for two dimensional lattices are clearly superior. 
In fact, the error bar on the exponents for two dimen- 
sional lattices is typically half that of one dimensional 
lattices. 

It is known that the Manna Model suffers from signif- 
icant logarithmic corrections [20] . Manna himself noted 
a "considerable curvature" in what should have been a 
straight line in a double logarithmic plot [11]. Similarly, 
Liibeck and Heger [56] found a surprising splitting of the 
Manna universality class in one dimension, which might 
also be due to the presence of significant corrections. 

The results for the finite size scaling exponents in Ta- 
ble IV suggest that one-dimensional systems are more dif- 
ficult to fit, with the alternative fitting functions Eq. (13) 
and Eq. (12) both clearly performing worse than in two 
dimensions. One might think that some of the prob- 



lems are caused by having much higher accuracy in the 
estimates in one dimension (and thus requiring more cor- 
rection term, as bad choices for the fitting function can 
no longer be hidden in a large statistical error), given 
that we spent, per lattice, typically about five times 
more CPU time than in two dimensions. The opposite is 
the case (probably because correlation times grow like a 
power law of the linear extent which are very large in one 
dimension) : Depending on the observable, relative errors 
are between a factor 2 and 10 worse in simple chain com- 
pared to square lattice. This holds similarly for other 
lattices, except for the nnn chain, on which we spent less 
than 1/6 of the CPU time we spent on the other one 
dimensional lattices. 

In general the relative statistical error vanishes like the 
inverse square of the CPU time, so that the product of 
the two gives a measure of "efficiency" of a lattice. Us- 
ing that measure, the simple chain is the most efficient, 
followed by nnn chain, rope ladder and finally the Futat- 
subishi lattice, which is about a factor 4 less efficient (4 
times the CPU time is needed for results with similar rel- 
ative error). This statement, however, is put in perspec- 
tive, by noting that we used variety of different hardware 
throughout. The two-dimensional lattices fall roughly in 
three classes: nc diagonal square lattice, Kagome lat- 
tice, jagged lattice, square lattice, followed by triangular 
lattice, Archimedes lattice and Mitsubishi lattice, and fi- 
nally the honeycomb lattice. The latter is clearly the 
worse (again by about a factor 4), while the triangular 
lattice is in the first group for some of the observablcs. 

There is a caveat, however. Within a given amount of 
CPU time and for a given system size, the nc diagonal 
square lattice (Fig. 2(d)) produces a larger number of 
avalanches, which are typically much smaller than those 
for other lattices, because the fraction of virtual neigh- 
bours is about 1.5 times higher for the nc diagonal square 
lattice than for, say, the honeycomb lattice, so that parti- 
cles are dissipated more frequently As a result, statistics 
are comparatively better for nc diagonal square lattice, 
which, however, may pose higher demands on correction 
terms with generally larger amplitudes. We could, how- 
ever, not identify a systematic behaviour in this respect. 
For example, the moments of jagged lattice are, within 
error, the same as for squarelattice, yet the latter has 
a much smaller average fraction of dissipative links. In 
fact, as discussed below, the same leading order ampli- 
tudes are found for the avalanche area distribution across 
all lattices of the same dimension. 

The one dimensional lattice perform particularly badly 
for the moment ratios. Essentially only those for the 
area size distribution can be fitted well and then pro- 
duce fairly consistent results (except for the nnn chain). 
Originally we expected improved scaling behaviour with 
the introduction of next nearest neighbour interaction, as 
it prevents degeneracy issues (and, say conserved quanti- 
ties) as they are sometimes observed on the simple chain 
[25, 66]. The two dimensional lattices, with the excep- 
tion of the triangular lattice, generally behave much bet- 
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TABLE VII: Comparison of the results in the present work to the estimates found in the literature. Many of the works quoted 
below have studied variants of the Manna Model. The values taken from [52, 53] and [54] in two dimensions are for the Oslo 
Model. The exponents marked as DP are those for the directed percolation universality class. They are derived via scaling 
laws [55]. 
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ter, when fitting moment ratios. Given the importance 
of the boundary conditions, it is remarkable how well 
all universal quantities addressed in the present work 
are reproduced, even when some two dimensional lat- 
tices like the Kagome lattice (Fig. 3(b)) and the jagged 
lattice (Fig. 2(b)) have rather complicated boundaries 
(although the latter is the square lattice in the thermo- 
dynamic limit) and many of the lattices have an aspect 
ratio of unity only asymptotically. 

To our surprise, the amplitudes A a ^ n in Eq. (15) ob- 
tained when fitting the moments of the area distribution 
seem to be universal themselves. According to Eq. (7) 



avalanche areas in one dimension follow 



A 



Jo 



dyy n - T *G a (y) 



(25) 



which is universal provided the metric factors a a and 
b a are, which is not normally the case. Since D a = d, 
however, the amplitude b a of the cutoff b a L Da is dimen- 
sionless. It is therefore reasonable to assume that it is 
universal. There is no reason, however, to assume that 
the same should hold for a a — the argument that a a is 
determined by normalisation does not hold as Eq. (5) ap- 
plies only asymptotically, i.e. the fraction of small event 
sizes which do not follow simple scaling can, in principle, 
vary from lattice to lattice. Moreover, a a is dimensionful 
since r ^ 1. Nonetheless, it turns out that it does not 
vary. As a result, to leading order, the moments of the 
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and in two dimensions 
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(32) 
(33) 
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as a function of the number N of sites, independent of 
the lattice type. We need to qualify the statement, by 
pointing out that in one dimension the amplitudes across 
different lattices are rather noisy, and in two dimension 
the amplitude of (a 1 ) has a goodness of fit of just under 
0.1 (we are using the [•] notation again here), while the 
other results in two dimensions all have a goodness of fit 
better than 0.5. 

Remarkably, had we fitted the area moments against a 
lattice-dependent multiple of N 1 / 2 , such as the perceived 
linear extent of the lattice, then this multiplier would 
have shown in the resulting amplitude A a ^ n . and so the 
apparently universal behaviour would not have come to 
light. 
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As explained above, the fitting function is a hypothe- 
sis and ultimately has an impact on the results. As the 
number of free parameters increases so does the suscepti- 
bility of the result on the initial condition. The approach 
described above, using fairly large systems (with weaker 
corrections) with not too small error bars, in conjunction 
with simple fitting functions, the initial values of which 
are determined by those with fewer terms, seems to pro- 
duce robust and reliable results. Comparing the results 
based on Eq. (15) and Eq. (13) in Table IV to those on 
the basis of Eq. (12) indicates that the former are supe- 
rior. An acceptable goodness of fit is reached for Eq. (12) 
only for those exponents that coincide within a bit more 
than one standard deviation with the estimates based 
on Eq. (15). Eq. (13), on the other hand, in summary 
(Table IV) coincides with Eq. (15), but some, individual 
finite size scaling exponents, such as [i[ s \ but also 



and D a , were estimated too poorly. 

In summary, the present work confirms the Abelian 
Manna Model as an SOC model that displays non-trivial, 
robust, reproducible, universal scaling behaviour in one 
and two dimension across different lattices. 
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